function [zz,rkhatsU,rkhatU] = find_gam_p(gam_p,RbU,LU,PU,lam_p,delta_p,sigrk_p)

rkhatsU = log(RbU*(1-1/LU)*(1+lam_p*(1-gam_p))-1+delta_p);

% Note: norminv(p) = (-sqrt(2))*erfcinv(2*p);

rkhatU = rkhatsU - sigrk_p*(-sqrt(2))*erfcinv(2*PU);
murkU = rkhatU;

% Note: int_{rkhatsU}^{infty} e^{rk}dF(rk) = exp(mu+sig^2/2)*normcdf((mu+sig^2-rkhatsU)/sig)
%       normcdf(x) = (1/2)*(1+erf(x/sqrt(2))).

pi = 3.141592653589793;
ErkU = exp(murkU+sigrk_p^2/2)*(1/2)*(1+erf((murkU+sigrk_p^2-rkhatsU)/(sigrk_p*sqrt(2))));
zz = RbU*(1-PU) - (ErkU+(1-delta_p)*(1-PU) - lam_p.*RbU^(2).*(1-gam_p).*(1+lam_p.*(1-gam_p)).*(LU-1)./(LU^2)./(RbU.*(1-1/LU).*(1+lam_p.*(1-gam_p))-1+delta_p).*(1./(sqrt(2.*pi).*sigrk_p)).*exp(-1/2.*((rkhatsU-murkU)./sigrk_p).^2));

end